#Replication Paper
#Gov 2001
#Anna Hopper & Sam Imlay
#March 27, 2013

#(ALSO SEE IMPUTING CODE FILE FOR SECOND HALF OF CODE)

setwd("~/Desktop")
setwd("~/Dropbox/Replication Paper")

#load data

#data<-read.csv("~/Desktop/PGoutsideTrim.csv")
data<-read.csv("~/Desktop/PGreplication.csv")

head(data)
#class(data$aq16_)
#install.packages("survival") 
library(survival)

# Subset and clean data: have to change the variable that lists 0-strongly dislike instead of just using integers 
newdata<-subset(data, touse==1)
head(newdata)
newdata$aq16_=as.character(newdata$aq16_)
newdata$aq16_=ifelse(newdata$aq16_=="", NA, newdata$aq16_)
newdata$aq16_=ifelse(newdata$aq16_=="0-strongly dislike", 0, newdata$aq16_)
newdata$aq16_=ifelse(newdata$aq16_=="10-strongly like", 10, newdata$aq16_)
newdata$aq16_=as.numeric(newdata$aq16_)
newdata$aq16_=ifelse(newdata$aq16_=="0-strongly agree", 0, newdata$aq16_)


####### Table 1 Model REPLICATIONS #######

# Table 1, Model 1: Distance 
#aq16_ is the party feeling thermometer, drdistp is physical distance from the constituency, and the conservative and liberal democrat variables are dummies)
clogit1.1<-clogit(votechoice~conservatives+liberal_democrats+aq16_+drdistp+strata(id)+cluster(aconsnam), method="approximate", data=newdata)
clogit1.1

#Table 1, Model 4 
clogit1.4<-clogit(votechoice~conservatives+liberal_democrats+aq16_+drdistp+ inccand+absdepdist+strata(id)+cluster(aconsnam), method="approximate", data=newdata)
clogit1.4

####### Simulations ########
#extract betas and sigmas for simulation from multivariate normal 
 
mean<-coef(clogit1.1)
vcv1.1<-vcov(clogit1.1)

mean4<-coef(clogit1.4)
vcv1.4<-vcov(clogit1.4)


###set values for covariates (this uses means for covariates and sets the values of distance)
LabRows<-subset(newdata, newdata$partyvote==1)

#Conservative matrix: 
ConsRows<-subset(newdata, newdata$partyvote==2)

#Liberal Democrat matrix: 
LDRows<-subset(newdata, newdata$partyvote==3)

##### Following Arzheimer and Evans:
##### Simulating from Model 4, using means and overall mean distance for the two other parties

#randomly draw betas 
library(MASS)
set.seed(02138)
betas<-mvrnorm(1000, mean4, vcv1.4) 

#### Run through for loops for each simulation in Table 3:

#### "Real" ####

#empty vectors for pi values to go in 
piL<-NULL
piC<-NULL
piLD<-NULL

#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
	# Insert real mean distances for each party
		j<-mean(LabRows$drdistp, na.rm=T)
		k<-mean(ConsRows$drdistp, na.rm=T)
		l<-mean(LDRows$drdistp, na.rm=T)
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))

	#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piL[i]<-catL/sum
		piC[i]<-catC/sum
		piLD[i]<-catLD/sum
}

mean(piL)
mean(piC)
mean(piLD)

mat<-cbind(piL, piC, piLD)
mean.point<-cbind(mean(piL), mean(piC), mean(piLD))
mat.2<-rbind(mat, mean.point)

# Coloring points by predicted probabilities of party vote
color<-matrix(nrow=nrow(mat), ncol=1)
for (i in 1:nrow(mat)){
	ifelse(max(mat[i,])==mat[i,1], color[i]<-"red", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,2], color[i]<-"blue", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,3], color[i]<-"yellow", color[i]<-color[i])
}
head(color)
color.2<-cbind(color, "black")


# Ternary graph!
ternaryplot(mat.2, scale=1, pch=16, cex=.3, col=color.2, dimnames=c("Labour", "Conservative", "Liberal Democrat"), dimnames_position="corner", main="'Real': Candidates at party distance means")
#grid_legend(0.8, 0.7, c(19, 19, 19, 19),
#	c("red", "blue", "yellow", "black"), c("Labour", "Conservative", "Liberal Democrat", "Mean"))
	
	

#### "Scenario 1" #### (distances held to mean(newdata$drdistp))

#empty vectors for pi values to go in 
piL<-NULL
piC<-NULL
piLD<-NULL

#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
	# Insert real mean distances for each party
		j<-mean(newdata$drdistp, na.rm=T)
		k<-mean(newdata$drdistp, na.rm=T)
		l<-mean(newdata$drdistp, na.rm=T)
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))

	#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piL[i]<-catL/sum
		piC[i]<-catC/sum
		piLD[i]<-catLD/sum
}

mean(piL)
mean(piC)
mean(piLD)

mat<-cbind(piL, piC, piLD)
mean.point<-cbind(mean(piL), mean(piC), mean(piLD))
mat.2<-rbind(mat, mean.point)

# Coloring points by predicted probabilities of party vote
color<-matrix(nrow=nrow(mat), ncol=1)
for (i in 1:nrow(mat)){
	ifelse(max(mat[i,])==mat[i,1], color[i]<-"red", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,2], color[i]<-"blue", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,3], color[i]<-"yellow", color[i]<-color[i])
}
head(color)
color.2<-cbind(color, "black")


# Ternary graph!
ternaryplot(mat.2, scale=1, pch=16, cex=.3, col=color.2, dimnames=c("Labour", "Conservative", "Liberal Democrat"), dimnames_position="corner", main="Scenario 1: Candidates at mean distance")



#### "Scenario 2" #### (Tories at 120km, others at mean(newdata$drdistp))

#empty vectors for pi values to go in 
piL<-NULL
piC<-NULL
piLD<-NULL

#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
	# Insert real mean distances for each party
		j<-mean(newdata$drdistp, na.rm=T)
		k<-120
		l<-mean(newdata$drdistp, na.rm=T)
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))

	#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piL[i]<-catL/sum
		piC[i]<-catC/sum
		piLD[i]<-catLD/sum
}

mean(piL)
mean(piC)
mean(piLD)

mat<-cbind(piL, piC, piLD)
mean.point<-cbind(mean(piL), mean(piC), mean(piLD))
mat.2<-rbind(mat, mean.point)

# Coloring points by predicted probabilities of party vote
color<-matrix(nrow=nrow(mat), ncol=1)
for (i in 1:nrow(mat)){
	ifelse(max(mat[i,])==mat[i,1], color[i]<-"red", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,2], color[i]<-"blue", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,3], color[i]<-"yellow", color[i]<-color[i])
}
head(color)
color.2<-cbind(color, "black")


# Ternary graph!
ternaryplot(mat.2, scale=1, pch=16, cex=.3, col=color.2, dimnames=c("Labour", "Conservative", "Liberal Democrat"), dimnames_position="corner", main="Scenario 2: Conservatives at 120 km")


#### "Scenario 2.1" #### (Tories at various distances, others at mean(newdata$drdistp))

dists<-seq(0, 120, by=5)

p.ests.cons<-matrix(data=NA, ncol=length(dists), nrow=1000)

for (m in 1:length(dists)){
# empty vectors for means and bounds

#empty vectors for pi values to go in 
piL<-NULL
piC<-NULL
piLD<-NULL

# Insert real mean distances for each party
		j<-mean(newdata$drdistp, na.rm=T)
		k<-dists[m]
		l<-mean(newdata$drdistp, na.rm=T)
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))


#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
		#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piL[i]<-catL/sum
		piC[i]<-catC/sum
		piLD[i]<-catLD/sum
		
	# Capturing mean Conservative vote share plus upper and lower 	bounds
	p.ests.cons[i,m]<-piC[i]
}

}


pdf("ConsVSD.pdf")
plot(dists, apply(p.ests.cons,2,mean), ylim=c(.05, .65), col="white", main="Conservative vote share by distance", xlab="Conservative candidate distance in km", ylab="Vote share")
lines(dists, apply(p.ests.cons, 2, mean), col="blue", lwd=2)
lines(dists, apply(p.ests.cons, 2, quantile, .025, col="blue"), lty=2)
lines(dists, apply(p.ests.cons, 2, quantile, .975, col="blue"), lty=2)
abline(h=.5, lty=2, col="red")
legend("topright", c("Mean", "95% confidence bounds", "50% probability"), lty=c(1, 2, 2), lwd=c(3, 2, 2), col=c("blue", "black", "red"))
dev.off()

#Without 50% line 
plot(dists, apply(p.ests.cons,2,mean), ylim=c(.05, .65), col="white", main="Conservative vote share by distance", xlab="Conservative candidate distance in km", ylab="Probability")
lines(dists, apply(p.ests.cons, 2, mean), col="blue", lwd=2)
lines(dists, apply(p.ests.cons, 2, quantile, .025, col="blue"), lty=2)
lines(dists, apply(p.ests.cons, 2, quantile, .975, col="blue"), lty=2)
#abline(h=.5, lty=2, col="red")
legend("topright", c("Mean", "95% confidence bounds"), lty=c(1, 2), lwd=c(3, 2), col=c("blue", "black"))
first.diffs<-p.ests.cons[,1]
min(first.diffs)
plot(density(first.diffs))
mean(first.diffs)
quantile(first.diffs, .025)
quantile(first.diffs, .975)


#### "Scenario 3" #### (LibDems at 120km, others at mean(newdata$drdistp))

#empty vectors for pi values to go in 
piL<-NULL
piC<-NULL
piLD<-NULL

#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
	# Insert real mean distances for each party
		j<-mean(newdata$drdistp, na.rm=T)
		k<-mean(newdata$drdistp, na.rm=T)
		l<-120
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))

	#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piL[i]<-catL/sum
		piC[i]<-catC/sum
		piLD[i]<-catLD/sum
}

mean(piL)
mean(piC)
mean(piLD)

mat<-cbind(piL, piC, piLD)
mean.point<-cbind(mean(piL), mean(piC), mean(piLD))
mat.2<-rbind(mat, mean.point)

# Coloring points by predicted probabilities of party vote
color<-matrix(nrow=nrow(mat), ncol=1)
for (i in 1:nrow(mat)){
	ifelse(max(mat[i,])==mat[i,1], color[i]<-"red", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,2], color[i]<-"blue", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,3], color[i]<-"yellow", color[i]<-color[i])
}
head(color)
color.2<-cbind(color, "black")


# Ternary graph!
ternaryplot(mat.2, scale=1, pch=16, cex=.3, col=color.2, dimnames=c("Labour", "Conservative", "Liberal Democrat"), dimnames_position="corner", main="Scenario 3: Liberal Democrats at 120 km")



#### "Scenario 3.1" #### (Lib Dems at various distances, others at mean(newdata$drdistp))

dists<-seq(0, 120, by=5)

p.ests.ld<-matrix(data=NA, ncol=length(dists), nrow=1000)

for (m in 1:length(dists)){
# empty vectors for means and bounds

#empty vectors for pi values to go in 
piL<-NULL
piC<-NULL
piLD<-NULL

# Insert real mean distances for each party
		j<-mean(newdata$drdistp, na.rm=T)
		k<-mean(newdata$drdistp, na.rm=T)
		l<-dists[m]
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))


#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
		#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piL[i]<-catL/sum
		piC[i]<-catC/sum
		piLD[i]<-catLD/sum
		
	# Capturing mean Conservative vote share plus upper and lower 	bounds
	p.ests.ld[i,m]<-piLD[i]
}

}


pdf("LDvsd.pdf")
plot(dists, apply(p.ests.ld,2,mean), ylim=c(.05, .65), col="white", main="Liberal Democrat vote share by distance", xlab="Lib-Dem candidate distance in km", ylab="Vote share")
lines(dists, apply(p.ests.ld, 2, mean), col="orange", lwd=2)
lines(dists, apply(p.ests.ld, 2, quantile, .025, col="black"), lty=2)
lines(dists, apply(p.ests.ld, 2, quantile, .975, col="black"), lty=2)
legend("topright", c("Mean", "95% confidence bounds"), lty=c(1, 2), lwd=c(3, 2), col=c("orange", "black"))
dev.off()


#### "Scenario 4" #### (Labor at 120km, others at mean(newdata$drdistp))

#empty vectors for pi values to go in 
piL<-NULL
piC<-NULL
piLD<-NULL

#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
	# Insert real mean distances for each party
		j<-120
		k<-mean(newdata$drdistp, na.rm=T)
		l<-mean(newdata$drdistp, na.rm=T)
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))

	#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piL[i]<-catL/sum
		piC[i]<-catC/sum
		piLD[i]<-catLD/sum
}

mean(piL)
mean(piC)
mean(piLD)

mat<-cbind(piL, piC, piLD)
mean.point<-cbind(mean(piL), mean(piC), mean(piLD))
mat.2<-rbind(mat, mean.point)

# Coloring points by predicted probabilities of party vote
color<-matrix(nrow=nrow(mat), ncol=1)
for (i in 1:nrow(mat)){
	ifelse(max(mat[i,])==mat[i,1], color[i]<-"red", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,2], color[i]<-"blue", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,3], color[i]<-"yellow", color[i]<-color[i])
}
head(color)
color.2<-cbind(color, "black")


# Ternary graph!
ternaryplot(mat.2, scale=1, pch=16, cex=.3, col=color.2, dimnames=c("Labour", "Conservative", "Liberal Democrat"), dimnames_position="corner", main="Scenario 4: Labour at 120 km")





#### "Scenario 4.1" #### (Labor at various distances, others at mean(newdata$drdistp))

dists<-seq(0, 120, by=5)

p.ests.lab<-matrix(data=NA, ncol=length(dists), nrow=1000)

for (m in 1:length(dists)){
# empty vectors for means and bounds

#empty vectors for pi values to go in 
piL<-NULL
piC<-NULL
piLD<-NULL

# Insert real mean distances for each party
		j<-dists[m]
		k<-mean(newdata$drdistp, na.rm=T)
		l<-mean(newdata$drdistp, na.rm=T)
	# Establish X values
xL<-c(mean(LabRows$conservatives, na.rm=T), mean(LabRows$liberal_democrats, na.rm=T), mean(LabRows$aq16_, na.rm=T), j, mean(LabRows$inccand, na.rm=T), mean(LabRows$absdepdist, na.rm=T))

xC<-c(mean(ConsRows$conservatives, na.rm=T), mean(ConsRows$liberal_democrats, na.rm=T), mean(ConsRows$aq16_, na.rm=T), k, mean(ConsRows$inccand, na.rm=T), mean(ConsRows$absdepdist, na.rm=T))

xLD<-c(mean(LDRows$conservatives, na.rm=T), mean(LDRows$liberal_democrats, na.rm=T), mean(LDRows$aq16_, na.rm=T), l, mean(LDRows$inccand, na.rm=T), mean(LDRows$absdepdist, na.rm=T))


#loop through simulated betas to get simulated 1000 simulated pi values for each party.

for (i in 1:1000){
		#numerator of pi 
		catL<-exp(xL%*%betas[i,])
		catC<-exp(xC%*%betas[i,])
		catLD<-exp(xLD%*%betas[i,])
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		piL[i]<-catL/sum
		piC[i]<-catC/sum
		piLD[i]<-catLD/sum
		
	# Capturing mean Conservative vote share plus upper and lower 	bounds
	p.ests.lab[i,m]<-piL[i]
}

}


pdf("LabVSD.pdf")
plot(dists, apply(p.ests.lab,2,mean), ylim=c(.05, .65), col="white", main="Labour vote share by distance", xlab="Labour candidate distance in km", ylab="Vote share")
lines(dists, apply(p.ests.lab, 2, mean), col="red", lwd=2)
lines(dists, apply(p.ests.lab, 2, quantile, .025, col="black"), lty=2)
lines(dists, apply(p.ests.lab, 2, quantile, .975, col="black"), lty=2)
legend("topright", c("Mean", "95% confidence bounds"), lty=c(1, 2), lwd=c(3, 2), col=c("red", "black"))
dev.off()






















#subset data to exclude places where all three candidates live less than 10 km 
newdata$close<-0

for (i in 1:887){
	k<-i*3
	ifelse(newdata$drdistp[k]<5 || newdata$drdistp[k-1]<5 || newdata$drdistp[k-2]<5, newdata$close[k]<-1, newdata$close[k]<-0)
	ifelse(newdata$drdistp[k]<5 || newdata$drdistp[k-1]<5 || newdata$drdistp[k-2]<5, newdata$close[k-1]<-1, newdata$close[k-1]<-0)
	ifelse(newdata$drdistp[k]<5 || newdata$drdistp[k-1]<5 || newdata$drdistp[k-2]<5, newdata$close[k-2]<-1, newdata$close[k-2]<-0)
}
table(newdata$close)

#new subsets
CloseData<-subset(newdata, newdata$close==1)
FarData<-subset(newdata, newdata$close==0)

#test new subsets to see if results change 
clogitC1<-clogit(votechoice~conservatives+liberal_democrats+aq16_+drdistp+strata(id)+cluster(aconsnam), method="approximate", data=CloseData)
clogitF1<-clogit(votechoice~conservatives+liberal_democrats+aq16_+drdistp+strata(id)+cluster(aconsnam), method="approximate", data=FarData)
clogitC1
clogitF1

#test for model dependence
clogit4<-model.frame(clogit1.4)
dim(clogit4)
dim(newdata)
#install.packages("WhatIf")
library(WhatIf)
counterfact4<-clogit4

for (i in 1:2078){
ifelse(counterfact4$liberal_democrats[i]==1, counterfact4$drdistp[i]<-120, counterfact4$drdistp[i]<-26)
}
counterfact4$conservatives<-mean(counterfact4$conservatives)
counterfact4$liberal_democrats<-mean(counterfact4$liberal_democrats)
counterfact4$aq16_<-mean(counterfact4$aq16_)
counterfact4$absdepdist<-mean(counterfact4$absdepdist)
counterfact4$inccand<-mean(counterfact4$inccand)
d<-whatif(~conservatives+liberal_democrats+aq16_+drdistp+inccand+absdepdist, data=clogit4, cfact=counterfact4) 



###########
# MISCELLANEOUS
###########
# Taking coefficients from Model 1 regression
mean.betas<-coef(clogit1.1)

# Extracting relevant x-values from data...
LabRows<-subset(newdata, newdata$partyvote==1)
xL<-(cbind(LabRows$conservatives, LabRows$liberal_democrats, LabRows$aq16_, LabRows$drdistp))
#Conservative matrix: 
ConsRows<-subset(newdata, newdata$partyvote==2)
xC<-(cbind(ConsRows$conservatives, ConsRows$liberal_democrats, ConsRows$aq16_, ConsRows$drdistp))
#Liberal Democrat matrix: 
LDRows<-subset(newdata, newdata$partyvote==3)
xLD<-(cbind(LDRows$conservatives, LDRows$liberal_democrats, LDRows$aq16_, LDRows$drdistp))

# Calculating probabilities
#numerator of pi 
		catL<-exp(xL%*%mean.betas)
		catC<-exp(xC%*%mean.betas)
		catLD<-exp(xLD%*%mean.betas)
	#denominator of pi
		sum<-catL+catC+catLD
	#pi values for each party 
		prob.L<-catL/sum
		prob.C<-catC/sum
		prob.LD<-catLD/sum

L<-(prob.L)
C<-(prob.C)
LD<-(prob.LD)

# Coloring by incumbent candidate or not...
length(newdata$inccand)
inccand<-matrix(nrow=887, ncol=1)
for (i in 1:887){
	k<-i*3
	j<-sum(newdata$inccand[k]+newdata$inccand[k-1]+newdata$inccand[k-2])
	inccand[i,1]<-j
}
head(inccand)

inccolor<-matrix(nrow=887, ncol=1)
for (i in 1:887){
	ifelse(inccand[i,1]==1, inccolor[i]<-"red", inccolor[i]<-"black")
}
inccolor

mat.0<-cbind(L, C, LD, inccand, inccolor)
dim(mat)
mat<-na.omit(mat.0)
mat.p<-na.omit(cbind(L, C, LD))

ternaryplot(mat.p, scale=1, pch=16, cex=.3, col=mat[,5], dimnames=c("Labour", "Conservative", "Liberal Democrat"), dimnames_position="corner", main="Predicted probabilities of party vote, Model 1")
grid_legend(0.7, 0.7, c(19, 19),
	c("red", "black"), c("Incumbent", "No incumbent"))
	

#######
# Coloring by inccand and incumbent's party
length(newdata$inccand)
inccolor<-matrix(nrow=887, ncol=1)
for (i in 1:887){
	k<-i*3
	inccolor[i,1]<-"gray"
	ifelse(newdata$inccand[k]==1, inccolor[i,1]<-"yellow", inccolor[i,1]<-inccolor[i,1])
	ifelse(newdata$inccand[k-1]==1, inccolor[i,1]<-"blue", inccolor[i,1]<-inccolor[i,1])
	ifelse(newdata$inccand[k-2]==1, inccolor[i,1]<-"red", inccolor[i,1]<-inccolor[i,1])
}
head(inccolor)


mat.0<-cbind(L, C, LD, inccand, inccolor)
dim(mat)
mat<-na.omit(mat.0)
mat.p<-na.omit(cbind(L, C, LD))

ternaryplot(mat.p, scale=1, pch=16, cex=.3, col=mat[,5], dimnames=c("Labour", "Conservative", "Liberal Democrat"), dimnames_position="corner", main="Predicted probabilities of party vote, Model 1")
grid_legend(0.8, 0.7, c(19, 19, 19, 19),
	c("gray","red", "blue", "yellow"), c("No Incumbent","Labour Incumbent", "Conservative Incumbent", "Liberal Democrat Incumbent"))



# Coloring points by predicted probabilities of party vote
color<-matrix(nrow=442, ncol=1)
for (i in 1:nrow(mat)){
	ifelse(max(mat[i,])==mat[i,1], color[i]<-"red", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,2], color[i]<-"blue", color[i]<-color[i])
	ifelse(max(mat[i,])==mat[i,3], color[i]<-"yellow", color[i]<-color[i])
}
head(color)


# Ternary graph!
install.packages("vcd")
library(vcd)
ternaryplot(mat.p, scale=1, pch=16, cex=.3, col=color, dimnames=c("Labour", "Conservative", "Liberal Democrat"), dimnames_position="corner", main="Predicted probabilities of party vote, Model 1")
grid_legend(0.8, 0.7, c(19, 19, 19),
	c("red", "blue", "yellow"), c("Labour", "Conservative", "Liberal Democrat"))
	


